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Abstract 

Sulci are localized furrows on the surface of soft materials that form by a compression-induced 
instability. We unfold this instability by breaking its natural scale and translation invariance, and 
compute a limiting bifurcation diagram for sulcfication showing that it is a scale-free, sub-critical 
nonlinear instability. In contrast with classical nucleation, sulcification is continuous, occurs in 
purely elastic continua and is structurally stable in the limit of vanishing surface energy. During 
loading, a sulcus nucleates at a point with an upper critical strain and an essential singularity in the 
linearized spectrum. On unloading, it quasi-statically shrinks to a point with a lower critical strain, 
explained by breaking of scale symmetry. At intermediate strains the system is linearly stable but 
nonlinear ly unstable with no energy barrier. Simple experiments confirm the existence of these two 
critical strains. 
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Sulci are usually seen in combination with complementary protrusions known as gyri 
on the surface of the primate brain, but are also seen on the palm of our hand, in our 
elbows and knees, in swollen cellular foams (such as bread) and gels, and in geological 
strata; a few representative examples are shown in Fig. [T] While observations of sulci 
are ancient, their systematic study is fairly recent; an early reference is the reticulation 
patterns in photographic gelatin [T], and there has been a small interest in these objects 
both experimentally j2HS] and theoretically [3j EH9], starting with the pioneering work of 
Biot j7] over the past 50 years. Despite this, there is no careful analysis of the fundamental 
instability and bifurcation that leads to sulci. Here we study the formation of a sulcus in a 
bent slab of soft elastomer, e.g. PDMS: as the slab is bent strongly, it pops while forming 
a sulcus that is visible in the lower right panel of Fig. [TJ releasing the bend causes the 
sulcus to vanish continuously, in sharp contrast with familiar hysteretic instabilities that 
pop in both directions. We find that sulcification is a fundamentally new kind of nonlinear 
subcritical surface instability with no scale and a strongly topological character, yet has no 
energetic barrier relative to an entire manifold of linearly stable solutions. We also argue 
that sulcification instabilities are relevant to the stability of soft interfaces generally, and 
provide one of the first physical examples of the consequences of violating the complementing 
condition [12] (during loading) and quasiconvexity at the boundary j6] (during unloading), 
keystones in the mathematical theory of elliptic partial differential equations and the calculus 
of variations. 

To understand the unusual nature of sulcification, we first recall Biot's calculation for the 
linear instability of the surface of an infinite half-space of an incompressible elastomer that 
is uniformly stressed laterally. Because the the free surface is the softest part of the system, 
and there is no characteristic length scale in the equations of elasticity or in the boundary 
conditions, instability first arises when the Rayleigh surface wave speed vanishes. For in- 
compressible rubber, Biot [7J showed that this threshold is reached when the compressive 
strain exceeds 45.3% |7J, at which value, all surface modes are unstable while the fastest 
growing one has an infinite wave number. Since every free surface looks like a half-space 
locally, Biot's instability lurks at every free boundary. Finite geometries typically break this 
infinite degeneracy and lead to a hierarchy of ordinary buckling instabilities that preempt 
the surface instability [TH]- However, since Biot's calculation was limited to a linear analysis 
of the problem, it could not address the question of whether the instability was supercritical 
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FIG. 1: Examples of sulci in (a) a primate brain (b) the arm of an infant, (c) sliced bread 
under lateral compression, and (d) the bent slab of PDMS used in our experiments. The results of 
a numerical simulation shown in red capture the form of the sulcus in the gel, as described in the 
text, with no adjustable parameters. The scale bar represents 2.3 cm in (a) 5 cm in (b) 2.5 cm in 

(c) and 0.33 cm in (d). 



or subcritical or its ultimate nonlinear saturation. Since the basic problem is scale free and 
translation invariant (the sulcification instability can arise anywhere along the surface), the 
nonlinear problem is numerically intractable without explicit regularization and a careful 
limiting process requiring that we unfold the sulcus literally and figuratively. 

Therefore we consider the bent strip geometry shown in Fig. [T]i, and break scale invari- 
ance by assuming that a thin skin of a stiff material is attached to the surface of the bent 
slab. Furthermore, the curvature maximum at the bottom of the horseshoe where the high- 
est strains are achieved naturally breaks translation invariance. For planar deformations, 
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the total energy of the system is given by 
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where x (X) is the deformation of a strip occupying a rectangular material volume Q C R 2 
and subject to the incompressibility constraint det (Vx) = 1, /x is the shear modulus of the 
incompressible elastomer, Bq is the stiffness of a semi-flexible skin and s is the arc length 
parameter of the upper surface T C dfl. Our model corresponds to having a simple neo- 
Hookean elastomer free energy for the bulk and a Bernoulli-Euler curvature energy for the 
skin. To understand the onset of the sulcification instability, we extremized the energy above 
using a custom-built finite element method with continuous strains and a hierarchical mesh 
(see Supplementary Information (SI)). We enforced incompressibility and self-contact using 
pressure fields and by assuming left-right symmetry about the sulcus. This model has three 
relevant length scales: a regularization length l r = \/ j^, the length of the strip L s , and its 



thickness W s . We use Hq and L s to scale all quantities so that B = and the aspect ratio 



of the strip L s /W s are the only dimensionless parameters. 

Our simulations start with an initially flat rubber strip that is bent and quasistatically 
compressed between parallel, rigid plates separated by a distance A. As the control pa- 
rameter A is decreased, compressive strain on the inner surface of the strip increases and 
ultimately drives sulcification. To ensure that the scale of the furrow is not numerically 
under-resolved, we use a recursively refined finite element mesh near the incipient sulcus to 
keep the mesh scale roughly an order of magnitude smaller than l r (See SI, Sec. B). Using 
a novel continuation method for variational inequalities (see SI, Sec. C), we computed both 
stable and unstable extrema of E (x), and explored the limit B — > 0. This yields the central 



result of our study, the family of bifurcation diagrams shown in Fig, 2(a) , where we plot the 
minimum height of the slab h as a function of A. 

Each h — A curve is a bifurcation diagram for a different value of B: solid lines represent 
linearly stable solutions, while the dotted lines represent linearly unstable solutions. Each 
curve follows the characteristic S-shape of a hysteretic transition, associated with a sudden 
change in h and formation or relaxation of a finite size sulcus when a critical value of A 
is passed in loading or unloading. For every value of B, extrema on the top branch have 
smooth surfaces, those on the middle branch have a pendant of size l r [see inset a of Fig. 
2(a)| , and those on the bottom branch have a self-contacting sulcus (insets b and c). As B 
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FIG. 2: (a) Bifurcation diagrams for the bent strip geometry shown in Fig. 1 (d), showing the 

scaled height h as a function of the scaled compression A, the bifurcation parameter for 

B £ [lCP 7 — 1CP 12 ] (yellow-magenta); solid lines mark stable equilibria and dotted lines mark 

unstable equilibria; thick gray lines highlight the asymptotic T-shaped diagram delimited by the 

Biot point on the left and the T point on the right. Solutions on the upper stable branches have 

smooth surfaces. Insets a-c show the structure (color=strain energy) for typical solutions on the 

other branches, rescaled to a fixed size, (b) The upper two curves show compressive strain in 

terms of the principal stretch at the eventual location of the sulcus X = on the smooth 

branch for values of A at the right and left fold points, respectively, for each B. The lowest curve 

shows the smallest value 1 — \ x attained on the surface outside the sulcus for the solution at the 

right fold point. The dashed lines correspond to the upper-critical strain predicted by Biot (upper 

line), and the new lower-critical strain extrapolation from our simulations (lower line), (c) The 

energy barrier AE for sulcification as a function of B sampling the solution at even intervals 
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between the left (bottom curve) and right fold (top curve) points. 



is decreased, (B G [1CT 7 — 10~ 12 ], l r /L s G [4.6 x 1CT 3 — 1CT 4 ]), the hysteresis in a typical 
loading cycle becomes atypically one-sided as the branch of unstable extrema (dotted lines) 
swings up toward the top stable branch and the S-shaped bifurcation diagrams converge to 
a master T-shaped diagram traced by the thick gray line with two critical points. 

For fixed A and B an unstable solution represents a saddle point in the energy landscape; 
its energy relative to the configuration with a flat surface is AE = E u (A) — E s (A) where 
E u (A) and E s (A) are the energies of the unstable and top, stable branches at A respectively, 
is an upper bound on the height of the barrier to nucleating a sulcus. Fig. 2(c)| shows AE 
as a function of B and confirms the convergence of the family of bifurcation diagrams 
toward the limiting T-shaped bifurcation diagram, as well as the existence of a nonlinear 
surface instability with no energetic barrier over an extended range of strains. We see 
that the instability thus differs significantly from traditional first order phase transitions 
in that the deformation is continuous, occurs in simple elastic continua and is well defined 
in the limit of vanishing surface energy. The presence of metastable region bracketed by 
a pair of critical strains along which the stable and unstable solutions coincide as the skin 
becomes vanishingly thin naturally explains the discrepancy between Biot's prediction and 
a large number of experiments on creasing and sulcification |3] (and references therein) 
over the past half century. Unfolding the instability without breaking translation symmetry 
at the surface, e.g. in a swollen, adhered layer of gel, then naturally leads to extreme 
sensitivity to imperfections, and a hierarchy of complex subcritical instabilities connecting 
Biot's instability and buckling (see SI, sec. E), and the ability to control sulcification [6J. 

As B — > 0, the sequence of saddle-node fold points encountered during loading converges 
to a limiting, infinitely sharp fold point when the surface strain at the lowest point on the 
inner surface of the horseshoe, X = 0, reaches a critical value of 45.6% consistent with 
Biot's classical result; in Fig. 2(b)| we see the convergence of the critical compressive stain, 
1 — \ x where \ x is the principal stretch of the deformation gradient Vx along the free 
surface, for finite B to Biot's predicted value at B — 0. A numerical linearized spectral 
analysis of the loaded slab also confirmed Biot's prediction that the Rayleigh surface wave 
speed vanishes at X = just as the critical strain is achieved, and corresponds to the 
failure of the complementing condition jH [12], wherein infinitesimal periodic solutions at 
the boundary grow at a rate that diverges as the inverse of the wavelength. Over-damped 
dynamical simulations-which trace steepest descent contours of the energy landscape-reveal 
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that nonlinear effects reorganize these surface waves into a self-similar furrowing process; 
after a short transient, depending on B, the growth of a sulcus, which occurs via rolling, 
not snapping, is described by the self-similar form x s (X, A*) + \^Xt\* (^K/^/xtj where A 
is a dimensional constant and v* is the numerically computed scale-invariant form of the 
sulcus [I], and x s (X, A) is the branch of smooth solutions and A* is the value of A at Biot's 
limiting fold point. 

The complementary sequence of saddle-node fold points for decreasing B encountered 
during unloading are actually "corners" associated with the loss/gain of a self-contacting 
sulcus as the unstable surface pendant just closes to form a cavity of fixed size l r . As 
these corners converge to the limiting "T-point", the maximal surface strain outside the 
self-contacting region approaches a limiting value of 35.4% that is attained at a sequence of 
points converging to X = 0. The convergence to this strain is traced by the lowest curve 



in Fig. 2(b) with the asymptote marked by the dashed line. The middle solid curve of Fig. 
2(b)| is another estimate of the critical strain computed by measuring the strain at X = 



for a sequence of extrema for corresponding values of A on the top branch. The T-point 
critical strain (like the Biot critical strain) is universal for free surfaces of incompressible 
materials, consistent with recent experimental observations [HE]; however they both change 
with applied normal stress (i.e. indentation), for compressible materials [4J etc. 

To understand why the T-point bifurcation and the entire unstable manifold are not 
captured by linearized analysis, we note that before the sulcus reaches the regularization 
scale l r , it shrinks according to the form x s (X, A) + I (A) (X// (A)) where / (A) > 



vanishes at the T-point and vy ~ v* (See insets b and c in Fig. 2(a) , and the relative scale 
factor of 6.5.) Since the elastic stress is determined by Vx, this transformation shrinks the 
size of the sulcus without altering the local stress balance; therefore all the material and 
contact non-linearities remain relevant even for vanishingly small sulci. 

We tested our theory with experiments using a commercial Sylgard 180 Elastomer to form 
36 x 26 x 4 mm slabs that were placed between parallel rigid plates attached to linear motors 
and compressed in small increments of 200 fim in a second, separated by 50s to allow for the 
equilibration of the slab. We tracked sulcification optically by imaging the refracted image 
of a laser sheet that passed through the slab along its bending axis. When the sulcus formed 
it sharply refocused the laser sheet into an almond-shaped caustic pattern surrounding a 
dark shadow (SI, sec. D). Fig. [3] (left) shows the evolution of a central raster scan of the 
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FIG. 3: A comparison of experimental (left) and computed (right) light intensity patterns cast 
by the central section of a laser sheet illuminating the slab's bending axis; the horizontal axis is 
transverse to the sheet and the vertical axis is the step number in a loading-unloading cycle, and 



corresponds to changes in A as shown in Fig. 2(a) 4 = formation of sulcus, when a caustic 



suddenly forms, 7 = vanishing of sulcus 

caustic pattern during a loading cycle (vertical axis) (see SI, sec. D). Analogous ray traced 
light distributions for the simulation, using the measured laser profile and assuming left-right 
symmetry, the measured system geometry, and B = 1CT 11 (physically l r ~ 18 /im), are also 
shown in Fig. fright) for comparison; the numbered red dots correspond to the numbered 
red dots in Fig. |2(a)[ We see that with no free parameters we can capture the one-way 
hysteretic transition associated with sulcification. 

The emergence of the T and Biot points, and intervening metastable region in the B — > 0, 
can be understood in terms of a nonlinear generalization of Biot's half-space problem. All 
our simulations show that when a half-space of incompressible elastomer is compressed by 
34.5% it has an infinite degeneracy of energy minimizers: the trivial flat configuration, and 
a continuous family of isolated sulci which are stable up to translation and rescaling (i.e. 
v (X) — > Iv (X//) for any number I > 0), i.e., these symmetries are spontaneously broken. 
Sulcification exchanges compressive strain for rotation and shear which are ultimately lim- 
ited by self-contact. Beyond the lower critical strain, forming a sulcus of size 0(1) <C 1 
releases energy over a region of size 0(l 2 ), equivalent to the failure of quasiconvexity at the 
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boundary. The spatial variations of Vx near the compressive strain maximum at X = 
act as symmetry-breaking perturbations and determine the ultimate scale I of the surface 
fold. When the compressive strain is not localized to a point, the size of the sulcus is not 
set by the local geometry of Vx and the domain, and the T bifurcation is sensitive to de- 
tails and potential interactions between multiple sulci resulting in reticulation [H |3], or in a 
combination of buckling and sulcification [T7| IT8] etc. (see SI, sec. E). 

More generally, T bifurcations might arise in elastic systems with internal interfaces and 
nucleat ionlike processes in elliptic systems where nonlinearities enter in a scale-free way, 
e.g. the formation of cavities, bubbles and cracks [HI [TJj. These processes are notoriously 
difficult to control, displaying extreme sensitivity to imperfections, and are associated with 
a discontinuous transition in the microscopic state characterized by a critical size nucleus; 
e.g, a bubble or crack will grow only once it has reached a threshold size. T points should 
exist in these systems in the limit when the surface energy vanishes and the size of the 
"defect" also vanishes, but the ratio of the two which corresponds to a critical pressure or 
stress remains finite. 

Acknowledgments: We thank D. Cuvelier for help with the experiments and the Kavli 
Institute and the Harvard MRSEC for support. 
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Supporting Information for "Unfolding the sulcus" 
by E. Hohlfeld and L. Mahadevan 



Here we provide more details about our simulation techniques, the full form of the in- 
equalities governing the self-contact problem inside the sulcus, our finite element discretiza- 
tion, the details of the rough continuation technique and experimental methods together 
with further experimental results, and a sketch of some generalizations to multiple buckling, 
sulcification and dynamics. 



In this first subsection we present, in detail, the equations we have used to model sulci- 
fication. The equations governing the interior of the rubber block, in dimensionless form, 
are 



energy density used in the main text; Aij := dxi/dXj. 

The contact problem for the surface within the fold of the sulcus is governed by a system 
of variational inequalities. We have assumed left-right symmetry for simplicity. Using a 
coordinate system where the 1-direction is parallel to the initially flat upper surface of the 
strip and the 2-direction is the outward normal to this surface, on the right half of the upper 
surface the inequalities describing contact are 



The first inequality expresses the restriction of the upper surface to the right half plane - 
by symmetry, this is the same as forbidding self-penetration. The second inequality is dual 
to the first and restricts the contact pressure, 7r, to be positive. The final equality is the 
Karush-Kuhn- Tucker condition which states that the contact pressure is zero if the upper 



Equations 




where S is the nominal stress tensor, p is the pressure, and W = \tr (A T A — 2) is the strain 



X\ > 0, 7T > 0, X\TX = 0. 
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surface is not in contact. A similar system of inequalities governs the interaction of bottom 
surface of the bent strip with the "apparatus" that bends the strip to drive sulcification on the 
upper surface. The regularizing bending stiffness term appears in the boundary conditions 
for the upper surface as 

The effects of the skin become important to the shape of the sulcus when the sulcus 
reaches the size ~ B3. Using the above formula and the scaling form of the sulcus at the 
T-point we see that 

1 _ 4 
S i2 ~ B x £?3 x B 3~ const. 

That is, a critical stress must emerge in the limit that the surface energy and "defect" size 
simultaneously go to zero as explained in the text, but the T-point stress/strain need not be 
the limiting stress computed here which is a function of position and the form of the surface 
energy. 



Finite element discretization. 



To simulate sulcification in the bent strip geometry, we discretized the equations of the 
incompressible neo-Hookean model using the finite element method. (See e.g. [8j for a 
detailed description the finite element method.) An example mesh for the finest scale sim- 
ulations presented in main text is shown in Supplementary Fig. |4| The elements in the 
mesh alternatively represent displacement degrees of freedom as tensor products of Her- 
mite polynomials or pressure degrees of freedom as quadratic Lagrange polynomials. This 
combination allows for continuous strains and stresses, which is important since the sulcifi- 
cation and Biot instabilities have characteristic strains; Hermite displacement elements are 
also required for compatibility with the surface bending energy. The contact pressure was 
described by continuous quadratic Lagrange surface elements. Because left-right symmetry 
was assumed, only the right half of the mesh shown was used, resulting in approximately 
9000 degrees of freedom for the finest scale simulation in the main text. 
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FIG. 4: The mesh used for the simulations in the main body of the paper is recursively refined 

on the surface near the sulcus. Displacement variables have continuous derivatives and are 
modeled as tensor product Hermite elements. Pressure variables are continuous and modeled as 

quadratic Lagrange elements. 



Rough continuation method 

We solved the discretized equations using rough continuation, which is a variant of pseudo- 
arclength continuation adapted for variational inequalities. Pseudo-arclength continuation 
is described in j3], it is Newton's method applied to a system with a control parameter in 
which the control parameter is treated as an additional variable on equal footing with the 
state variables. Before explaining rough continuation, we will first abstractly describe the 
usual pseudo-arclength continuation method. 

Let the state of some system be described by the vector x G M", let A be a control 
parameter, and suppose the system is governed by the equation 

f(x,A) = 0, (1) 

where f G lR n . At fixed A, a Newton's method iteration is defined as 

x fc+1 = x A - Df (xfc, A) -1 (f (x fe , A)) 

where Df is the n x n derivative matrix of f . Pseudo-arclength continuation computes the 
connected component of solutions to (JlJ) whose closure includes some initial point (xo, Ao) 
rather that just a single solution for some specific value of A. An iteration in this method is 

Xfc+i = x,fc — PfcLfc (f (xfc, Afc)) 
Afc+i = Afc + 5k 
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where Lk is any left inverse of the n x (n + 1) matrix 
Pk is the projection operator 



<9f 

L>f (x fc , A fe ) , — (x fc , A fc ) 



in which t k £ IR n+1 is the unit row vector orthogonal to all the rows A fe , and 8 k is the 
A-component of t&. At the end of N iterations at step r, we make a new initial guess for 
(xI+W 1 ) by 

{x[ +1 ,\ r l +1 ) = (x r N ,\ r N ) + At r N 

where A is a small, usually adaptively chosen step size. Each iteration of the method 
accomplishes an approximate projection back to the solution curve x(A). The method has 
an ambiguity in the sign of A. This is fixed by setting t r N ■ t^ +1 > 0. Pseudo-arclength 
continuation inherits the quadratic rate of convergence of Newton's method whenever the 
matrix A (x, A) is continuous in x and A. 

Both Newton's method and pseudo-arclength continuation can easily handle constraints 
posed as Lagrange multipliers, such as the constant volume constraint of our incompress- 
ible strip. Inequality constraints, like contact, pose a problem, however. Linear inequality 
constraints can be posed as 

(C • (x - R)) . > 0, (p).>0, (C • (x - R)). (p). = 

for a fixed vector R and matrix C where p is the vector of Lagrange multipliers, and 
the notation (p) i means the i th component. The active set method can be used to solve 
these inequalities at fixed A. In this method, a set of "active" constraints are guessed and 
the inequality constraints in the active set are replaced with equality constraints while the 
remaining inequality constraints are ignored. A Newton iteration is then attempted. If 
the result of the Newton iteration violates other constraints, these are turned on for the 
next iteration, while if a Lagrange multiplier becomes negative, the associated constraint is 
turned off in the next iteration. 

Rough continuation is a synthesis of the pseudo-arclength continuation method with the 
active set method that solves two interrelated problems. First, whenever a constraint turns 
on or off, the matrix A must change discontinuously. And second, the sign of A is again 
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ambiguous since each new Lagrange multiplier changes the sign of det Df . In practice, these 
are only problems at parameter values where a constraint is just turning on or off along the 
length of the curve x(A), which results in a discretization induced corner in the otherwise 
smooth curve x (A). We can solve both problems by simply replacing Newton iterations with 
active set iterations in the pseudo-arclength method except when the set of active constraints 
is not converging even while the iterations remain bounded. In this case we attempt to guess 
the correct "forward" direction on the underlying smooth curve (i.e. without discretization), 
and then add a small multiple of tj~ to (xfc,Afc) at each iteration so to "nudge" (xjt+i,Afc + i) 
in this direction. Hence the resulting projection step is not always orthogonal to the current 
tangent t^, but rather orthogonal to our best guess for the chord connecting (x^y, X r N ) and 
(x^ +1 ,AJ +1 ). This modification drives the iterations away from where a constraint is just 
turning on or off so that the active set is fixed and ordinary pseudo-arclength iterations can 
resume. In practice the method easily handles large steps in which several constraints turn 
on or off at once. 



Experiments 

This subsection presents details of the experimental data presented in the main text 
and additional results from other experiments that help to characterize sulcification. These 
additional results include observations of the caustic pattern cast by a fully three dimensional 
sulcus, observations of the dynamic formation of a sulcus, and measurements of the force 
applied to the confining plates that bend the PDMS slab into a horseshoe shape and drive 
sulcification. 

As stated in the main text, we compared the predicted shape of the sulcus to the ex- 
perimental shape by making detailed a comparison between refracted light patterns from 



laser illumination. Supplementary Fig. 5(a) shows ray tracing for a bent slab at maxi 



mum compression. The green lines passing through the cylindrically curved inner surface of 
the horseshoe are simply refracted as if by a cylindrical lens, while the green lines passing 
through the sulcus are refracted into a dominant caustic pattern surrounding a dark shadow. 
Light passing through the ends of the horseshoe is refracted in to a secondary pattern of 
caustics (blue lines) and also into a set of whispering gallery modes (red lines). 

A corresponding light pattern for a three dimensional sulcus in the experimental geometry 
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is shown in Supplementary Fig. 5(b)| The sample is illuminated by a laser sheet along its 
symmetry axis (vertical axis in the figure). The unrefracted light which missed the sample 
can be seen as the bright blue bands in top and bottom extremities of the figure. The red line 
in the figure is the raster along which data was sampled at intervals in the compression cycle 
to make the Fig. 4 (left) of the main text. The outermost bright blue ring in Supplementary 



Fig. 5(b) corresponds to the green caustic pattern in Supplementary Fig. 5(a) , the left half 
of which was partially obscured by the experimental apparatus. The inner most bright blue 



ring in Supplementary Fig. 5(b) corresponds to the blue caustic pattern in Supplementary 



Fig. 5(a) The ring-like shape of the caustic patterns is a consequence of edge and three 
dimensional effects which cannot be described by our two dimensional model. 



Supplementary Fig. 5(c) shows representative observations of the dark shadow cast by a 
sulcus in another sample during a loading cycle. In frame (1), the cylindrical curvature of the 
bent PDMS strip broadened the laser sheet passing through the sample; the unrefracted sheet 
can be seen as the bright band at the top of each frame, and the dark gap immediately below 
this band is an artifact of the sample's edge. When the system achieved the critical strain 
for sulcification, a system spanning sulcus appeared and cast the almond shape shadow seen 
in frame (2). The shape of the shadow reflects that the sulcus tapered at its upper and lower 
extrema. Frame (3) shows the shadow pattern shrank in width and depth during unloading 
and suggests that the sulcus retracted along its length while simultaneously becoming more 
shallow. Frame (4) shows the initial configuration for the second compression cycle. The 
clearly visible dark line running down the center of the image suggests that self-adhesion 
may have prevented the sulcus from completely unfolding; this line is also partially visible in 
frame (3). We observed that the line gradually filled in over a period hours to days. When 
the line was present in the refracted light pattern, we also saw a faint line in the surface of 
an unfolded sample when we observed the sample at glancing incidence. 

Supplementary Fig. [6] shows how sulcification disturbs a layer of small glass beads dusted 
on the surface of another sample. The rolling motion of the forming sulcus plows the beads 
out of the self-contacting region, causing them to pile up in the opening of the furrow. When 
the sample is unfolded, the depth of the sulcus is made apparent by the width of the depleted 
region. We see that the depth of the sulcus changes gradually along its length, tending to 
zero at the sample's edges. This is consistent with the idea that the local strain state sets 
the (local) scale of the sulcus. 
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Using high speed video, we observed that the sulcus and corresponding caustic pattern 
nucleated at point near the center of bent strip and then grew to span the width of the bent 
strip. We tried to quantify this nucleation process by analyzing high speed movies of the 
refracted light pattern, but the detailed patterns and growth speeds were highly variable. 



A typical result is shown in Supplementary Fig. 5(d) , in which the concentric circular lines 
trace the outline of an observed dark shadow in the refracted light field at 20 ms intervals; 
color also indicates elapsed time in seconds. The vertical axis is along the length of the 
sulcus. 



Supplementary Figs. 7(a) and 7(b) show equilibrium and time series measurements of the 
force applied to the pair of plates laterally confining the experimental sample. In Supple- 



mentary Fig. 7(a) , force measurements traced the lower branches of the green lines during 
compression until a sulcus formed and the force jumped to a less compressive value; the 
upper green curves were traced during unloading and we did not observe a second jump ex- 
pected for typical hysteretic processes. Black and red lines are computed from simulations 
for B = 1CT 11 , or l r pa 18 frni. We fitted the vertical scale and offset for the red curve to the 



corresponding section of the green curves. Supplementary Fig. 7(b) shows the time series of 
force data for the jumps in the green curves of Supplementary Fig. |7(a)| on successive com- 
pression cycles. Downward motion is an increment to A, upward motion is a sulcification 
event. We sampled the force reading at 100 Hz and then averaged in one second windows to 
produce the curves shown in the figure; the error bars are the standard deviation of the force 
in an averaging window. There is no apparent correlation between the plate motion (change 
in A) and the precise moment of sulcification. This is consistent with the local nature of the 
sulcification instability. However, each successive event occurs at lower applied load; this 
may result from imperfect unfolding of the sulcus due to self- adhesion. 



Generalizations: multiple creases, buckling and dynamics 

We have argued that nucleation processes, especially in hyperelastic solids with non- 
convex strain energies, can be understood via the existence of a scale-symmetry breaking 
solution to a scale-free auxiliary problem which is a nonlinear generalization of Biot's linear 
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shows an example of the ray tracing used to compute Fig. 4(right) of the main 
text. In Fig. 5(b)| a sulcus at maximum compression refracts a laser sheet into a caustic pattern; 
the bright bands at the top and bottom are the unrefracted laser sheet that missed the sample. 
The red line is the raster along which data were taken make Fig. 4(left) of the main text. Fig. 



5(c) shows refracted light from an another sample at four points in a compression cycle: (1) the 
initial configuration, (2) the first appearance of a sulcus, (3) as the sulcus quasistatically shrinks 
to a point, (4) the initial configuration after one cycle. The residual scar seen in (3) and (4) is also 



visible by eye at glancing angle and vanishes over a period of hours to days. Fig. 5(d) shows the 
evolving shape of the dark shadow cast by developing sulcus at 20 ms intervals (computed by 
thresholding frames from a high speed movie); color indicates time in seconds and the vertical 
axis is along the sulcus. In contrast with the quasi-static light distributions, a detailed three 
dimensional simulation is needed to calibrate thgse data, but we see that the sulcus rapidly grows 

from a point in a few tenths of a second. 



FIG. 6: Small glass beads dusted on another sample are scattered by the formation of sulcus 
revealing the self-contacting region when the sample is unfolded as the dark scar crossing the 
image. The depth of the sulcus changes gradually, tending to zero near the sample's edges, 
consistent with the idea that the local strain state controls the scale of the sulcus. 




experiment 
theory, compression 
theory, rarefaction 



12 13 14 15 16 17 18 
A (mm) 



19 20 21 22 




(a) 



FIG. 7: The green lines in Fig. 
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7(a) 



show the measured force-displacement curves for multiple 
loading cycles in the experiment discussed in the main text. The red and black lines are computed 

from the simulation with the vertical offset and scale of the compression branch (red) as fitting 
parameters. Fig. |7(b)| shows time traces of the force applied by the compressing plates during the 
compressive phase of the loading cycles: downward steps coincide with increments to A, upward 

steps indicate sulcification events. 

problem. In the context of nonlinear elasticity, this problem is to minimize the energy 



Ea,y (v) = I W (A + Vv) - W (A) cOC 

JFyr 



over all functions v e W l,p (i<Y,IR d ) where A is identified with the value of Vx at some 
point Yefl, Fy is all of space if Y is an interior point and Fy is a half-space if Y is a point 
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on the (smooth) boundary of Q, p is an exponent related to the growth of W at large A, 
and d is the dimension of space. Because of the scale-free nature of this auxiliary problem, 
i.e. E A Y is a homogenous function of degree d in L under the mapping 

v (X) -»■ Lv (X/L) , L > 0, 

a minimizer can only exist if 

/ W (A + Vv) - W (A) rfX > 0, Vv G C] (B r (x) , M d ) , r > 0, X G F Y . 

This latter condition embodies both quasiconvexity (when dist ^X, 8Fy J > r) and quasi- 
convexity at the boundary (otherwise). Quasiconvexity was invented by Morrey and is a 
necessary condition for the functional -Ea,y to be sequentially weakly lower semicontinuous 
on the Sobolev space W 1,p (-Fa,y,K^) , it is also sufficient for this property when W is con- 
tinuous and has polynomial growth in A [1, 2J. Sequential weak lower semicontinuity can be 
used to prove a minimizer exists when Ea,y is bounded from below. Quasiconvexity at the 
boundary was invented by Ball and Marsden as a necessary condition for a local minimizer in 
the topology of C 1 (-Fay? to be a local minimizer in the comparatively weaker topology 
of C° (F AtY ,R d ) Here it is needed to ensure that E A y is bounded from below. 

The relationship between quasiconvexity and quasiconvexity at the boundary for our 
auxiliary problem is analogous to the relationship between ellipticity (in the interior) and 
satisfaction of the Complementing Condition (i.e. ellipticity at the boundary) for quadratic 
energy functionals. These conditions result from linearizing the two quasiconvexity condi- 
tions at v = 0. 

The existence of a non-trivial minimizer (or minimizing sequence) for £a ,y determines the 
critical deformation gradient A for nonlinear instability in the original problem and implies 
that one of the quasiconvexity inequalities is saturated. When a non-zero minimizer exists, 
say v, there is a continuous family of such minimizers (related by scaling and translation), 
each of which solves the Euler-Lagrange equation, and so picking one is the spontaneous 
breaking of scale-symmetry in the auxiliary problem. Generally we then expect, as happens 
in our simulation, that E a ,y (v) = determines a hyper surface in the space of deformation 
gradients A G M. dxd and that -Ea,y (v) changes sign as this hypersurface is crossed. This 
implies that some form of quasiconvexity is just lost at the critical deformation gradient and 
corresponds to the onset of metastability in the original problem because the system develops 
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an instability toward nucleation (i.e. the scaling motion of v) while remaining linearly stable. 
However, we do not know if the nucleation process is itself sub-critical or super-critical with 
respect to the violation of quasiconvexity, nor can we anticipate the nucleation dynamics. 

When the variation of strain (and the domain shape at a boundary point) act to stabilize 
the nucleation instability, a T-bifurcation results and this bifurcating branch acts as a center 
manifold, guiding the evolution of the instability. Because these local variations act to keep 
the growing nucleus compact, the bifurcation and dynamics will be insensitive to the geom- 
etry at large distances from X . This is not necessarily true if the variations in strain and 
geometry are not stabilizing. In this case the nucleation instability becomes sub-critical with 
respect to the violation of quasiconvexity, and we can expect that the nucleation dynamics 
will become quite complicated. 

As a simple example, consider a uniformly compressed beam of incompressible rubber 
such as in the inset to Supplementary Fig. [8j this figure sketches a hypothetical bifurcation 
diagram for this system. In this free beam geometry, surface waves on opposite sides of 
the beam interact to break Biot's surface instability into a hierarchy of sub-critical buckling 
modes; a qualitative bifurcation diagram for the first buckling mode is traced by the blue 
line in the figure. Red crosses mark subsequent linear instabilities in the straight beam 
and tend to an accumulation point marking the Biot point. When the beam is compressed 
beyond the buckling threshold, a small initial disturbance grows exponentially in time. 

Quasiconvexity is still lost at the critical strain of ~ 35%, but could happen before or 
after the onset of buckling, depending on the aspect ratio of the beam. In this geometry, 
sulci nucleated at different points can interact with each other and with surface waves 
in a dynamical setting. In particular, the inset to Supplementary Fig. [8] suggests that 
buckling and sulcification can interact to enhance one another: sulcification lowers the beams 
effective bending stiffness while bending increases the local compressive strain which drives 
sulcification. The resulting sub-critical instability might be a kind of kinked buckling, but 
the precise location of the kinks, and so the precise development of the instability will be 
sensitive to imperfections and the initial conditions. In the simplest case when the instability 
is initiated by a localized perturbation near the midpoint of the beam and the strain in the 
initial state is above ~ 35%, we expect the instability initially resembles a nucleation process 
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governed by the scaling form 

x (X) +i(i - io )v( IIi | 7o) 

with the function L (t) determined by the dynamical law and the time to determined by 
initial conditions. For example if 

W = - V -^A (VX) ' 
which describes gradient descent, then L (t) = A 7 t 7 with 7 = 1/2. If 

<9 2 x „ dW ,„ , 

= V-T^r(Vx), 



dt 2 dA 

which describes inertial dynamics, then L (t) has the same form with 7 = 1. If the dynamics 
are viscoelastic, but local in time, e.g. if the dynamical law has the form 

v.s(vx,v|)^o, 

where S is the nominal strain, then L (t) = e*/ A . The green curve traces the bifurcation 
diagram for this process. Since kinked buckling is sub-critical, below the critical strain where 
quasiconvexity is lost, a small size sulcus will shrink according to the scaling dynamical law: 

Xo( x) + i(to - t )v( II 2L_). 

Notice that when the dynamics are either inertial or gradient descent, the sulcus will 
nucleate (or vanish) at a finite time, but that for viscoelastic material, the sulcus cannot form 
spontaneously (or completely vanish). Assuming that PDMS is governed by a viscoelastic 
dynamics, this potentially explains why the metastable portion of the computed bifurcation 
diagram is experimentally accessible and the existence of the vertical black line in SI Fig. 
2(c)(4). 
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FIG. 8: Here we show a schematic bifurcation diagram for the finite beam shown in the inset. 

This geometry destabilizes both the Biot-point and the T-point leading to a hierarchy of 
sub-critical buckling instabilities accumulating at the Biot-point (red crosses and blue line), and 
sub-critical sulcification (green line). Black arrow lines schematically show the dynamics as 

explained in the text. 
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